Introducción

El presente documento evidencia el proceso de obtención de las víctimas documentadas, imputadas y estimación del subregistro (desagregadas por tiempo) usando las funciones confirm_files, read_replicates, filter_standard_cev, summary_observed, combine_replicates, mse y combine_estimates del paquete verdata. Específicamente se replicará la figura (sup-izq) de la página 12 del informe metodológico.

¿Qué debemos tener en cuenta antes de usar el paquete?

Antes de mostrar este proceso, es importante definir algunos conceptos para que usted, como usuario del paquete, pueda no solo replicar la estructura del documento metodológico, sino también pueda explorar los datos e implementar nuevos análisis relacionados con el conflicto armado en Colombia.

Así, es importante tener en cuenta los siguientes conceptos:

A. Datos observados

Son aquellos datos de la víctima que ya estaban presentes en la base integrada, es decir, la base que fue consolidada a partir de la información de más de 100 bases de datos de diversas instituciones en Colombia. Sin embargo, si bien algunas víctimas cuentan con la información completa en variables como sexo, edad, etc., otras víctimas presentan campos faltantes en dichas variables, por ende, no tenemos certeza de cuáles fueron las características de las víctimas, es decir, tenemos y siempre tendremos incertidumbre sobre el valor o características reales de dichas personas.

B. Campos imputados

¿Cómo se aborda este tipo de vacío en la documentación? Se hace a través de un modelo estadístico de imputación múltiple que permite tener una base de datos donde ahora se observan valores posibles respecto a la edad de la víctima, presunto responsable, etnia, etc. A estas observaciones la denominamos campos o datos imputados (estadísticamente hablando. No confundir con imputación en términos jurídicos).

C. Estimación del subregistro

Si bien a través de la base integrada contamos con una cantidad considerable de víctimas registradas, es importante entender que aún faltan víctimas por documentar: aquellas que nunca fueron registradas en ninguna de datos. Para abordar el subregistro o los datos faltantes se utiliza una metodología que en estadística se llama Estimación por Sistemas Múltiples (ESM).

Autenticando e importando la base de datos (réplicas)

Se comienza autenticando e importando la base de datos de desaparición, esto a través de dos funciones de nuestro paquete: la función confirm_files y read_replicates. La primera función autentica los archivos que han sido descargados de cualquier sitio en internet en el que se encuentren. Considerando que cada violación tiene hasta 100 réplicas, esta función permite autenticar cada uno de estos archivos sin necesidad de leerlos a R. Esto, en caso de querer ahorrar recursos computacionales, o en caso de que no vaya a realizar su análisis con todas las réplicas. Esta función devolvera el mensaje “You have the right file!” si los archivos son iguales a los publicados, o el error “This file is not identical to the one published. This means the results of the analysis may potentially be inconsistent.” si no lo son.

verdata::confirm_files(here::here("verdata-replicas/verdata-parquet"), 
                                     "desaparicion", 1, 10)

Además, la función read_replicates permite 2 cosas: leer las réplicas a R en una sola tabla (ya sea a partir de un formato csv o parquet) y la segunda es verificar que el contenido de las réplicas sea exactamente igual al publicado. Cuando el argumento crash tiene su valor por default (TRUE), la función retorna un objeto (base de datos) si el contenido es igual, y el mensaje “The content of the files is not identical to the ones published. This means the results of the analysis may potentially be inconsistent.” si el contenido de la base fue previamente alterada/modificada, lo que quiere decir que los análisis que el usuario quiera hacer serán corruptos y conllevarán a resultados erróneos. Si usted quiere ver los datos en R, aunque el contenido haya sido alterado, puede cambiar el argumento crash a FALSE, y, en ese caso, podrá ver los datos, junto con el mismo mensaje de advertencia.

replicas_datos <- verdata::read_replicates(here::here("verdata-replicas/verdata-parquet"), 
                                     "desaparicion", 1, 10)

paged_table(replicas_datos, options = list(rows.print = 10, cols.print = 5))

Vemos pues que tenemos 1 720 970 registros, nuestras réplicas (o versiones de los datos) van desde la número 1 hasta la 10, vemos además que nuestros datos tienen información sobre la categoría de edad de la víctima, el presunto perpetrador, el sexo, el año del hecho, la pertenencia étnica, entre otros. Sin embargo, para centrarnos en un análisis más específico, tal como el realizado para el informe metodológico, procederemos a transformar y/o filtrar algunas variables.

Filtrando las réplicas acorde con el filtro del informe metodológico

La función filter_standard_cev nos permite transformar o filtrar nuestra información. Por ejemplo, aquellos víctimas que se documentaron como víctimas de la ex-guerrilla FARC-EP posterior a 2016, pasaron a ser víctimas de otras guerrillas, ya que este primer grupo oficialmente dejó de existir después de dicho año (perp_change = TRUE)

replicas_filtradas <- verdata::filter_standard_cev(replicas_datos,
                                                   "desaparicion", 
                                                   perp_change = TRUE)

paged_table(replicas_filtradas, options = list(rows.print = 10, cols.print = 5))

Víctimas documentadas

Después de aplicada la anterior función es momento de obtener la información de las víctimas documentadas por año del hecho. Estos datos en particular son aquellos que ya venian en la base integrada y que en ocasiones contenían campos faltantes en algunas de las variables. Usaremos pues la función summary_observed para calcular dicha documentación.

Como se observa, los argumentos de la función son: 1) la violación a analizar desaparicion; 2) nuestros datos replicas_filtradas; 3) strata_vars que en nuestro caso será la variable de yy_hecho; 4) le sigue el argumento de conflict_filter que filtra a aquellas personas que fueron víctimas dentro del marco del conflicto armado (variable is_conflict == TRUE) o no (variable is_conflict == FALSE).

Esta función también incluye un argumento denominado 5) forced_dis_filter el cual aplica únicamente a esta violación. Esta indica si la víctima fue desaparecida de forma “forzada” (forced_dis == TRUE) o no (forced_dis == FALSE). Para otras violaciones este argumento siempre será “FALSE”.

También contamos con otros argumentos: 6) edad_minors_filter que filtra por víctimas menores de edad (= TRUE) documentadas por los proyectos y/o instituciones; 7) include_props que permite incluir el cálculo de las proporciones para nuestras variables de interés (= TRUE); 8) prop_obs_na que permite incluir la proporción de NA’s en nuestra tabla (= TRUE) y 9) digits que es un argumento opcional en el cual podemos establecer el número de dígitos para redondear nuestros resultados (que por defecto es 2).

tabla_documentada <- verdata::summary_observed("desaparicion",
                                               replicas_filtradas, 
                                               strata_vars = "yy_hecho",
                                               conflict_filter = TRUE,
                                               forced_dis_filter = TRUE,
                                               edad_minors_filter = FALSE,
                                               include_props = FALSE,
                                               prop_obs_na = FALSE)

paged_table(tabla_documentada, options = list(rows.print = 4, cols.print = 4))

tabla_documentada <- tabla_documentada %>%
    mutate(yy_hecho = as.numeric(yy_hecho)) %>% 
    arrange(desc(observed))

g <- tabla_documentada %>%
    ggplot(aes(x = yy_hecho)) +
    geom_line(aes(y = observed, color = "Observado"),  size = 1) +
    theme_minimal() +
    theme(axis.text.x = element_text(size = 11, angle = 90),
          axis.title.y = element_text(size = 11),
          axis.ticks.x = element_line(size = 0.1)) +
    scale_x_continuous(breaks = seq(1980, 2016, 2)) +
    theme(legend.position = "bottom") +
    labs(x = "Año",
         y = "Número de víctimas",
         color = "") +
    scale_colour_manual(values = c("Observado" = "#434343"))

g

Posterior a esto se procede a aplicar la función de combine_replicates que, como su nombre lo indica, permite combinar las réplicas para obtener lo que denominamos “la media de la estimación puntual” junto con el intervalo de confianza que permite dimensionar la incertidumbre de la imputación. Para esta función se siguieron las reglas de combinación de Rubin, que, si desea ver con más detalle de qué se trata, el libro Flexible Imputation of Missing Data de Stef van Buuren aborda paso a paso este proceso.

Ahora, esta función se compone de los siguientes argumentos: la violación a analizar desaparicion; 2) tabla_documentada, es decir, el data frame derivado de la función summary_observed; 3) nuestra base de datos filtrada replicas_filtradas; 4) strata_vars que será nuevamente nuestra variable de año del hecho; 5) conflict_filter que filtra a aquellas personas que fueron víctimas dentro del marco del conflicto armado (variable is_conflict == TRUE) o no (variable is_conflict == FALSE).

Esta función también incluye un argumento denominado 5) forced_dis_filter el cual aplica únicamente a la violación de desaparición. Esta indica si la víctima fue desaparecida de forma “forzada” (forced_dis == TRUE) o no (forced_dis == FALSE). Para otras violaciones este argumento siempre será “FALSE”.

También contamos con otros argumentos: 6) edad_minors_filter que filtra por víctimas menores de edad (= TRUE); 7) include_props que permite incluir el cálculo de las proporciones para nuestras variables de interés (= TRUE); y 8) digits que es un argumento opcional en el cual podemos establecer el número de dígitos para redondear nuestros resultados (que por defecto es 2).

tabla_combinada <- verdata::combine_replicates("desaparicion",
                                                tabla_documentada,
                                                replicas_filtradas, 
                                                strata_vars = "yy_hecho",
                                                conflict_filter = TRUE,
                                                forced_dis_filter = TRUE,
                                                edad_minors_filter = FALSE,
                                                include_props = FALSE)

paged_table(tabla_combinada, options = list(rows.print = 10, cols.print = 5))

tabla_combinada <- tabla_combinada %>% 
     arrange(desc(imp_mean))

g <- tabla_combinada %>%
    mutate(yy_hecho = as.numeric(yy_hecho)) %>% 
    ggplot(aes(x = yy_hecho)) +
    geom_line(aes(y = observed, color = "Observado"),  size = 1) +
    geom_line(aes(y = imp_mean,  color = "Imputado"), size = 1) +
    theme_minimal() +
    theme(axis.text.x = element_text(size = 11, angle = 90),
        axis.title.y = element_text(size = 11),
        axis.ticks.x = element_line(size = 0.1)) +
    scale_x_continuous(breaks = seq(1980, 2016, 2)) +
    theme(legend.position = "bottom") +
    labs(x = "Año",
         y = "Número de víctimas",
         color = "") +
    scale_colour_manual(values = c("Imputado" = "#1F74B1", 
                               "Observado" = "#434343"))

g

Como se indicó, la línea de color negro muestra los datos documentados; esta documentación se refiere a víctimas de las cuales sabemos que efectivamente fueron víctimas dentro del marco del conflicto y que fueron desaparecidas de forma forzada; pero, ¿qué hay de las víctimas sin información acerca de estas características? Pues bien, esta información (junto con la ya documentada) se encuentra en la línea azul, la cual muestra las víctimas después del proceso de imputación múltiple. Es decir, como se indicó anteriormente, acá estamos incluyendo aquellas víctimas para las que no teniamos conocimiento de dicha información.

Así, después del proceso de imputación estadística y con un con un nivel de confianza del 95% se evidencia que hubo entre 9 250 a 9 336víctimas de desaparición forzada en el 2 002 con un promedio de 9 293. Es decir, esto implica que este promedio es la mejor estimación puntual respecto al número de víctimas de dicho año; no obstante, hay que tener en cuenta que siempre tendremos la incertidumbre de la imputación y que dicho fenómeno estará representado por el intervalo de confianza del 95%.

Procedemos a estratificar para el posterior proceso de estimación del subregistro de víctimas.

Proceso estratificación para estimaciones

Con el fin de controlar la heterogeneidad en las probabilidades de captura (ver más de este concepto en el informe metodológico procederemos a estratificar nuestra información (por año del hecho y la variable de “es desaparición forzada”)

replicas_estratos <- replicas_filtradas %>% 
  dplyr::mutate(is_conflict = as.integer(is_conflict)) %>% 
  dplyr::filter(is_conflict == 1) %>% 
  dplyr::filter(is_forced_dis == 1)

paged_table(replicas_estratos, options = list(rows.print = 10, cols.print = 5))

Seguido de esto procedemos a estratificar. Es importante que usted como usuario vea que este proceso es netamente artesanal, es decir, usted puede usar su propio código o funciones para realizar este proceso que, en nuestro caso, será a través de una función previamente creada (fuera del paquete verdata) para facilitar este ejercicio:

stratify <- function(replicate_data, schema) {
    
    schema_list <- unlist(str_split(schema, pattern = ","))
    
    grouped_data <- replicate_data %>%
        group_by(!!!syms(schema_list))
    
    stratification_vars <- grouped_data %>%
        group_keys() %>%
        group_by_all() %>%
        group_split()
    
    split_data <- grouped_data %>%
        group_split(.keep = FALSE)
    
    return(list(strata_data = split_data,
                stratification_vars = stratification_vars))

}

Entonces, en primera instancia creamos una función que necesita de dos (2) argumentos:

• La primera está relacionada con los datos que en nuestro caso son los replicas_estratos, es decir, la información del primer paso.

• La segunda son nuestras variables (o los nombres) de estratificación (schema). Todas estas variables se encuentran -precisamente- en nuestros replicas_estratos.

En términos generales, lo que hace esta función es: primero agrupa por nuestras variables de estratificación y guarda esta información en un objeto llamado grouped_data; también definimos el nombre de nuestra estratificación, es decir por nuestras variables de interés. Luego, split_data va a dividir nuestros datos previamente agrupados (en group_data), retornando una lista: strata_data que será lossubconjuntos de datos de acuerdo con la estratificación y stratification_vars que contiene las combinaciones (nombres) de nuestras variables.

replicas <- paste0("R", 1:10)
results <- list()

for (replica in replicas) {
  strata_data <- replicas_estratos[replicas_estratos$replica == replica, ]
  
  var <- NULL
  strata <- stratify(strata_data, c("yy_hecho", var, sep = ","))
  
  results[[replica]] <- strata
}

replica6_grupo6 <- results[["R6"]][["strata_data"]][[6]]

paged_table(replica6_grupo6, options = list(rows.print = 10, cols.print = 5))

replica6_grupo6_estrato <- results[["R6"]][["stratification_vars"]][[6]]

paged_table(replica6_grupo6_estrato, options = list(rows.print = 10, cols.print = 5))

Estimación víctimas por año del hecho

Ahora, definido nuestros estratos procedemos a calcular las estimaciones para esta violación en particular; acá entra nuestra función mse del paquete verdata. Esta función lo que nos permite es -por ejemplo- preparar los datos: como se indicó al principio, esta función tomará como insumo la información de las fuentes, es decir, aquellas columnas que comienzan por in_. También filtra por las fuentes válidas, es decir, si un estrato no es estimable nos arrojará NA lo que nos indica que ese estrato no cuenta con al menos 3 fuentes con al menos un (1) registro de violencia. Además -entre otros- nos permite revisar si tenemos estimaciones precalculadas, lo que nos permitirá ahorrar tiempo ya que estas provienen de procesos de estratificación que se realizaron en la CEV.

mse(
  stratum_data,
  stratum_name,
  estimates_dir = NULL,
  min_n = 1,
  K = NULL,
  buffer_size = 10000,
  sampler_thinning = 1000,
  seed = 19481210,
  burnin = 10000,
  n_samples = 10000,
  posterior_thinning = 500
)

Algunos de los argumentos de esta función se explican de la siguiente forma¹

• stratum_data: Data frame que incluye la información del estrato de interés (dataframes que creamos antes).

• stratum_name: Es el nombre del estrato.

• estimates_dir: Es la ruta (opcional) o el path de la carpeta o archivo que contiene las estimaciones precalculadas.

estimaciones_dir <- here::here("estimaciones")

Al final el resultado será un data frame con cinco (5) columnas: una columna que indica si el estrato es válido o no (TRUE o FALSE); el número de muestras N de la distribución posterior (NA si el estrato no es válido); las fuentes válidas que se usaron en la estimación valid_sources; el número de observaciones de las listas que son válidas del estrato (n_obs) y el nombre del estrato stratum_name. Si un estrato es estimable, este dataframe va a contener 1000 muestras para cada estrato. Esto lo veremos en un momento.

estimates <- list() 

replicas <- paste0("R", 1:10)

for (replica in replicas) {
    
    strata_data <- results[[replica]]$strata_data
    stratification_vars <- results[[replica]]$stratification_vars
    
    estimate_mse <- purrr::map2_dfr(.x = strata_data,
                                    .y = stratification_vars,
                                    .f = ~ {
                                      mse <- verdata::mse(stratum_data = .x,
                                                          stratum_name = .y,
                                                          estimates_dir = estimaciones_dir)
                               mse$replica <- replica
                               mse
                           })
    
    estimates[[replica]] <- estimate_mse
}

union_mse_tab <- do.call(rbind, estimates)

paged_table(union_mse_tab, options = list(rows.print = 10, cols.print = 8))

union_mse_tab <- arrow::read_parquet(here::here("dane-examples/output-estimacion/yy_hecho.parquet"))

union_mse <- union_mse_tab %>%
    mutate(stratum_name = paste(pull(stratum_name, 1),
                                sep = "-"))

paged_table(union_mse, options = list(rows.print = 10, cols.print = 8))

tabla_final <- union_mse %>%
    rename(replicate = replica) %>% 
    select(-validated, -valid_sources)

paged_table(tabla_final, options = list(rows.print = 10, cols.print = 5))

final_agrupacion <- tabla_final %>%
    group_by(stratum_name)

estimates_tabla <- final_agrupacion %>%
    group_split() %>%
    map_dfr(.f = combine_estimates) %>%
    bind_cols(group_keys(final_agrupacion)) %>% 
    select(stratum_name, N_025, N_mean, N_975) %>% 
    rename(yy_hecho = stratum_name)

paged_table(estimates_tabla, options = list(rows.print = 10, cols.print = 5))

Vemos pues que, después de agrupar, tenemos los resultados derivados de estimates_tabla cuyo código permite separar nuestras categorías (o o años), luego utilizamos la función de map_dfr en el que podemos aplicar combine_estimates a cada una de nuestras categorías. Por último -con bind_cols- unimos estos resultados con las variables de agrupación, es decir, con yy_hecho que la denominamos stratum_name.

En tal sentido, el proyecto “CEV-HRDAG-JEP” estimó que para el año 2002 el número de víctimas desaparecidas de forma forzada estuvo entre y 26 530 con una probabilidad del 95% (intervalo de credibilidad). En otras palabras, hay una alta probabilidad (95% de credibilidad) de que el número de de victimas en este año se encuentre dentro de este rango, mostrando a su vez el valor más probable de 21 850 victimas.

Por último uniremos esta base con nuestros resultados combinados, cortaremos la varianza de la estimación y graficaremos nuestros resultados:

estimates_final <- estimates_tabla %>% 
    mutate(yy_hecho = as.character(yy_hecho))

tabla_final <- dplyr::left_join(estimates_tabla, tabla_combinada)
    
paged_table(tabla_final, options = list(rows.print = 10, cols.print = 5))

n_warn <- 19000

combine_estimates_year <- tabla_final %>%
    mutate(max_var = case_when(
      (N_975 > n_warn) ~ n_warn,
      TRUE ~ NA_real_)) %>%
    mutate(N_975 = ifelse(!is.na(max_var), max_var, N_975))

combine_estimates_year <- combine_estimates_year %>% 
    arrange(desc(N_mean)) %>% 
    mutate(yy_hecho = as.numeric(yy_hecho))

mr_observed_ttl <- glue("Observado")
mr_replicates_ttl <- glue("Imputado")
mr_universo_ttl <- glue("Estimado")

g <- combine_estimates_year %>%
  ggplot(aes(x = yy_hecho)) +
  geom_line(aes(y = observed,
                fill = mr_observed_ttl), color = "black", size = 1) +
  geom_line(aes(y = imp_mean,  fill = mr_replicates_ttl),color = "#1F74B1", 
            show.legend = FALSE, size = 1) +
  geom_ribbon(aes(ymin = N_025, ymax = N_975, fill = mr_universo_ttl),
              alpha = 0.5) +
  geom_point(aes(y = max_var), pch = 21, color = 'darkgreen', fill = "green",
             size = 1, stroke = 2) +
  theme_minimal() +
  xlab("") +
  ylab("Número de víctimas") +
  theme(axis.text.x = element_text(size = 11, angle = 90),
        axis.title.y = element_text(size = 11),
        axis.ticks.x = element_line(size = 0.1)) +
  scale_x_continuous(breaks = seq(1985, 2020, 5)) +
  theme(legend.position = "bottom") +
  scale_fill_manual(values = c("darkgreen", "#1F74B1", "black" ), name = "")

print(g)